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Abstract The full set of published radial velocity data (52 measurements from Keck + 58 
ones from ELODIE + 17 ones from CORALIE) for the star HD37124 is analysed. Two 
families of dynamically stable high-eccentricity orbital solutions for the planetary system 
are found. In the first one, the outer planets c and d are trapped in the 2/1 mean-motion 
resonance. The second family of solutions corresponds to the 5/2 mean-motion resonance 
between these planets. In both families, the planets are locked in (or close to) an apsidal 
corotation resonance. In the case of the 2/1 MMR, it is an asymmetric apsidal corotation 
(with the difference between the longitudes of periastra 4 a) ~ 60°), whereas in the case of 
the 5/2 MMR it is a symmetric antialigned one {Aw = 180°). 

It remains also possible that the two outer planets are not trapped in an orbital resonance. 
Then their orbital eccentricities should be relatively small (less than, say, 0. 15) and the ratio 
of their orbital periods is unlikely to exceed 2.3 — 2.5. 

Keywords planetary systems ■ resonance ■ stability ■ periodic orbits • statistical methods 



1 Introduction 

For now, the star HD 37124 is believed to host three Jovian planets. The innermost planet 
'b' was discovered bv lVogt etal.l ( l2000l) . This planet moves on a low-eccentric orbit with a 
pe riod of Pt » 150 da ys. Soon after this, the s econd planet 'c' was discovered independently 
bv'Udry et al.l ( l2003h and lSutler et al. lllool). At that time, its mass and orbital parameters 
(e.g. period Pc ~ 2000 days) were highly uncertain. Finally, \^_gt et al. ( 2005) announced 
discovery of the third planet. Its orbit was most likely located between the orbits of planets 
'b' and 'c'Q 
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' There is no clear consensus between researchers about notation of planets in the system. The innermost 
planet is always denoted as 'b', but the notation for the outer pair of planets may vary. We use the same 
notation as used in The Extrasolar Planets Encyclopaedia by J. Schneider, www.exoplanet.eu. Namely, 
we denote the innermost, the outermost and the intermediate planet by the letters 'b','c','d', respectively. 
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Still, the radial velocity (RV) data are insufficient to obtain reliable estimations of pa- 
rameters of this system directly. Any attempt to obtain a best-fitting RV solution inevitably 
leads to dynamically unstable orbital configuration disintegrating after a very short time 
due to high eccen tricities of t wo ou ter planets. To force the fitti ng algorithm to find a sta- 
ble co nfiguration, IVogt et all ll2005h fixed the value of Cc at 0.2. iGozdziewski et al] ( I2006L 
l2008h presented a detailed analysis of the Keck RV data involving constraints of dynamical 
stability. The stable orbital configurations from these works span a very wide region, with 
ratio of orbital periods of the outer planets, Pc/Pd, ranging from ~ 2.1 to ~ 2.9. We have 
to ascertain that, in fact, the Keck RV data only outline a wide region of acceptable orbital 
configurations, whereas significant constraints on the orbits of these planets are set mainly 
by the stability requirement. 

The aim of the present paper is to describe the set of most likely orbital solutions for 
the system of HD37124, based on the analysis of the complete set of RV data published 
(incorporating Keck, ELODIE, and CORALIE measurements). The structure of the paper 
is as follows. In Section |2l the RV datasets used in the paper are described. In Section [3] 
the statistical methods used in the paper are discussed. In Sections |4l [5] and |6l the results 
of the analysis are presented. In Section [T] the dynamical behaviour of the resulting orbital 
configurations is considered. In Section [H the hypothesis of existence of an extra planet in 
the system is tested. 



2 Radial velocity data 

The most precise publicly available radial velocity data for HD37124 were published in the 
paper by Vogt et al. ( 2005). These 52 measurements were obtained at the Keck telescope and 
span about 8.4 yr (between 1996.9 and 2005.3). They show RV uncertainties from 2.1 m/s 
to 3.7 m/s. Also, the observations at ELODIE and CORALIE instruments were made. Al- 
though these are not so precise, they could significantly increase the temporal coverage of 
the full RV time series. Unfortunately, these data were not pub lished in a table f orm and only 
a graph of these m easurements is available in the paper by lUdrv et alj ( l2003i) . We apply a 
similar approach as lFerraz-Mello et alj ( l2005h and lSeauge et alj ( |2008|) used for HD82943. 



We recon struct the rneasure ments, their uncertainties, and their dates from the graph pub- 
lished in ( IUdrvetal.L[2003l) . The radial velocities themselves and their error bars can be 
reconstructed quite accurately (better than 1 m/s accuracy). The reconstructed Julian dates 
of observations have typical errors of ~ 1 day. This is admissible also, because the shortest 
orbital period Pi, w 150 days in HD37124 is much longer. The 58 reconstructed ELODIE 
data points cover about 7.2 yr between 1995.0 and 2002.2 and possess RV uncertainties 
from 7 m/s to 19 m/s. The 17 reconstructed CORALIE data points cover about 1.4 yr be- 
tween 1999.8 and 2001.2 and possess RV uncertainties from 6 m/s to 20 m/s. Thus the span 
of the combined time series is about 10.3 yr. This combined time series incorporate 7 = 3 
independent time series of Nj{i = 1,2,3) RV measurements Vji(i = 1,2, .. . ,Nj) having the 
'stated' RV uncertainties (Tmeas.^! and made at the epochs f^,. These datasets are plotted in 

Fig. [T] top panel. 

It was demonstrated in ( lBaluevil2008bl l3) that high-precision RV measurements in planet 
search surveys often suffer from periodic (annual) systematic errors which may originate 
from various sources. Sometimes, these systematic errors may reach the magnitude ~ 10 m/s, 
especially for data published several years ago, when data reduction algorithms have not 
been debugged to a perfect state yet. It is necessary to account for these systematic errors in 
our analysis. For this purpose, a simple harmonic model of these errors Acos(2;r(f — T)/lyr) 
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Fig. 1 Top panel: RV 
measurements for the star 
HD37124. Different RV 
offsets were assigned to the 
ELODIE, CORALIE, and 
Keck datasets, in order to 
separate them from each 
other. Bottom panel: RV 
curve for the be st stable fit 
from cGozdziewski et all 
120081) , its RV residuals, 
and its RV residuals after 
correction of the best-fitting 
sinusoidal annual drift in the 
ELODIE data. 



is adopted below. Here, the semi-amplitude A and the time shift t are the extra free param- 
eters to be determined from the time series. In the next section we will see that ELODIE 
data always show a significant annual drift of radial velocity with large semi-amplitude 
A ~ 20 m/s. Unfortunately, the number of CORALIE measurements is too small for reli- 
able modelling of their possible annual errors. Seemingly, it is better to avoid this modelling 
for the CORALIE data. The existence of significant annual errors in Keck measurements is 
uncertain. We will consider different models of the RV curve below, with and without the 
annual term in the Keck data. 

One might ask that since the data from ELODIE and CORALIE are less accurate, suffer 
from significant systematic errors, and increase the total time coverage by only 20— 25%, 
why to include them into analysis at all? If the orbits of planets in the system were deter- 
mined well by the Keck data alone then indeed addition of such RV data would not signifi- 
cantly improve the precision of the estimations. However, the Keck data alone provide only 
quite mild constraints on the orbits of the two outer planet. The major problems are due to 
ill-determinacy of the orbital period and eccentricity of the outermost planet. In this case, we 
should compare RV datasets in the sense of absolute rather than relative increase of the total 
time base. About ten ELODIE data points, which span two years before the regular obser- 
vations of HD37124 started at the Keck observatory, cover almost full orbital period of the 
intermediate planet and about one third or almost half of the orbital period of the outermost 
planet. During this time segment, the corresponding RV variations are about ~ 30 m/s for the 
planet d and ~ 10 — 20 m/s for the planet c. Such arcs of the RV curves are very important for 
constraining the corresponding orbital periods, which result in more strict constraining of 
the whole set of parameters. Of course, the Keck data remain the main source of information 
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and drive the fit, but the data from ELODIE are also important, because they can help to rule 
out a large fraction of inacceptable fits. This is i llustrated by the bottom p anel in Fig.[T] We 
can see that the orbital solution from the work dOozdziewski et all |2008|) fits satisfactorily 
all available RV data in the range 1997-2005, except for the ELODIE measurements made 
in 1995-1996, which show systematic deviation reaching 50 m/s. This makes the mentioned 
orbital solution significantly less credible. Other solutions, not being ruled out by the Keck 
data, may produce even larger deviations and thus could be easily ruled out even by very 
inaccurate measurements. However, data from different observatories have different statis- 
tical properties and it is necessary to merge them extremely carefully, in order to set correct 
statistical weights to different datasets. This problem will be considered in the next section. 

3 Statistical analysis: principles and definitions 

3.1 The general RV model and the system of parameters 

Let us write down the model of the radial velocity measurements obtained at 7* observatory 
at time t as 

i"j(f,P) =Mobsj(?,Pobsj)+M*(f,P*), j=l,2,...J. (1) 

Here, the full vector p of free parameters to be estimated consists of elements of vectors 
Pobsj(i = 1,2, and p*. The function fiobsj in O represent an observatory-specific 

part of the measured radial velocity: 

'J 

Mobs,>(f,Pobsj-) = CQj + A j„ cos{2n{t - Tj„)/P/„) (2) 

The constant velocity term co.j and parameters Aj„ , Tj„ of possible systematic errors form the 
vectors Pobs ,j of unknowns. The quantities Pj„ — the periods of the systematic errors — are 
assumed a priori known. In this paper, we will consider only the cases Sj = (no systematic 
errors) and Sj = 1 with Pj\ being the annual period. The function /i^. in Ul is the common 
radial velocity term incorporating RV signals due to unseen companions orbiting the star: 

r jY 

/i*(f,p*) = ^ c„f" + K„{cos{a„ + U,,) +e„cos(B„). (3) 

«— 1 17 = 1 

The coefficients c„ describe possible long-term polynomial (of degree r in general, we will 
consider only the cases r = 0, no trend, and r = 1, linear trend, below) trend in the RV data. 
This trend may be induced by possible distant unseen companions in the system with periods 
longer than the time span of the observations. Other terms in ^ represent Keplerian veloc- 
ities induced by planets {.yV = 3 in our case). The coefficients c„, RV semi-amplitudes 
K„ and the orbital elements A„ (the mean longitude at certain fixed epoch), e„ (the eccen- 
tricity), a>„ (the argument of the periastron), P„ (the orbital period) form the vector p* of 
planetary parameters to be estimated. The quantity i)„ in ((Sjl is the true anomaly of the «* 
planet (evidently, it depends on the time and on the parameters A„,P„,e„). 

The model ^ does not take into account interactions between the planets. This is ad- 
missible for the case of HD37124, because after ~ 10 yr of RV observations, the outermost 
planets 'c' and 'd' (which, as we will see below, represent the main source of dynamical ac- 
tivity in the system) have completed only two and four revolutions around the star On such 
a time scale, the planetary perturbations could be only significant in the case of sufficiently 
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Fig. 2 The difference between the Newto- 
nian N-body and multi-Keplerian RV model s 
for the fit from iGozdziewski et all 120081) . 
The clear growth of this difference up to 
~ 50 cm/s reflects the fact that the orbital el- 
ements were given for the epoch of the first 
Keck observation. If the reference epoch was 
closer to the middle of the observation win- 
dow, the RV difference would be bounded by 
20 — 30 cm/s over the whole time segment. 



close approaches between planets on high-eccentricity orbits. However, such configurations 
are unrealistic, because close approaches usually represent the source of instability and lead 
the system to disintegration. For example, the difference between the best-fitting Keplerian 
and N-body RV models (see Fig. [2| ) do no t exceed 50 cm/s for the best stable orbital con- 



figuration from dGozdziewski et alll2008h . For stable orbital solutions that are considered 
below, these deviations do not exceed ~ 30 cm/s and produce only ~ 10^"' relative error in 
the RV residuals r.m.s. Bearing this in mind, only the Keplerian RV models are used in the 
paper. They allow much faster computations than the N-body ones. 

The minimum mass of the planet m sin ; (here i is the orbital inclination to the sky plane) 
and its semi-major axis a can be derived as 

m sin i^K{ \ = ^KP^'I^mI'^ , 



\2TtG J 



-siP^I^Ml'^ , (4) 



where ^ w4.919- IQ-^ \Mjup-Mr:^l^ -m'^ -day^^/^ ■ s] and w 1.957-10-2 [AU-M^j^^^ • 
day-2/^] are constant factors, G is the gravitational constant. Mi, is the mass of the star, and 
K = K\/ 1 — e^ is the modified semi-amplitude. As it is well known, the inclination cannot 
be constrained using the Keplerian RV model. 

The errors of the approximate equalities in Q are about m/M* ~ lO^-^ and are much less 
than the statistical uncertainties of estimations given below. For example, the shift in the RV 
semi-amplitude of about m/M* ~ 10^^^ would lead to an error in the RV of only 1—3 cm/s 
and would be completely invisible in the resulting fit. Regular motions are not sensitive to 
such small shifts in the planetary masses. If the planets are not trapped in a mean-motion 
resonance, the same proposition holds true for the semi-major axes as well. However, we 
deal below with resonant motions, which are sensitive to the shifts of the semi-major axes of 
&(m/Mi,). For these cases, we need to provide more decimal digits (excessive with respect 
to the statistical uncertainties) for the values of the semi-major axes, in order to allow the 
reader to reproduce the results of long-term numerical integrations discussed below. For 
this purpose we should state the coord inat e system in which w e refer our estimations. It 
was noted bv lLissauer & Riveral ( l200ll) and lLee & Peald ( |2003|) that it is better to interpret 
the orbital parameters of the Keplerian model as osculating ones referenced in the Jacobi 
coordinate system. In the Jacobi coordinates, the osculating orbital period and semi-major 
axis of an fe* planet are connected by the relation 

ak = J^Pl'^ Ul, + Y,mj\ . (5) 



6 



Here we need to substitute the values of the planetary masses themselves, and thus to assume 
some values for orbital inclinations. As it will be discussed in Section[5] we assume that the 
planetary system is seen edge-on, that is = 90° . 



3.2 Estimations of the parameters 



To analyse the RV data described in Section |2l we need to estimate the so-called RV jitter 
a? which increases the full RV uncertainties as ai. 



full 



CT^ + Cmeas softens the differ- 



ences between the weights of observations cx: 1 / aj,,, . In the a strophysical part, this RV jitter 
is inspired by various activity in the st ar (e.g. 



Wrightll2005h . but often incorporates instru- 



mental effects as well. As is shown bv lBaluev ( 2008ch , the effective RV jitter may be quite 
different for different instruments, even for one and the same star. Therefore, we should 
perform the merging of RV datasets from different observatories very carefully, with correct 
assignment of statistical weights to these data. T o do it, we use here the maximum-likelihood 
approach described in the paper ( lBaluevl . [2008Q) . This algorithm includes a built-in estima- 
tion of the effective RV jitter (simultaneous with the estimation of usual parameters), which 
allows us not to rely on a low-precision astrophysical estimations of a^,. Moreover, this al- 
gorithm allows to perform a separate estimation of the effective RV jitters for the datasets 
from different observatories. This algorithm uses the maximization of the modified log- 
likelihood function of the A' RV observations vy,- (their errors are assumed to be uncorrelated 
and Gaussian-distributed) which is defined as 



In^: 



J 

EE 

>=i /=i 



r— 2 + In CTfuU,;,- 



-N\nV2n max. 



(6) 



full.j/ 



j + CTj^eas ji ij = 1, 2, ... ,7; I = 1,2,.. .Nj) and the correction divisor y = 



Here, a^ 

1 — d/N with d being the number of degrees of freedom in our RV model (i.e., the number 
of free parameters, d = dimp). The divisor y allows to perform a 'preventive' reduction of 
the statistical bias in the estimations of the RV jitter. Since we aim to use the new objective 
function instead of the usually used one, we need to intro duce a new mea sure of the 
goodness-of-fit, which would be based on In accordance with lBaluevI ( l2008d) . to assess 
the quality of a given orbital fit, the following goodness-of-fit measure is used below: 



e"" VV2;i: w 0.2420if 



-\/N 



(7) 



This function is measured in the same units as radial velocity (i.e., in m/s). It characterises 
the overall scatter of RV measurements around the model. However, to allow a comparison 
with previous works on HD37124, the traditional r.m.s. goodness-of-fit measure is used 
below as well. 



3.3 Assessing the reliability of orbital fits 

It is not enough to find an RV fit with a small scatter of res iduals. To interpret t he resulting 
estimations, we need to assess their reliability. It is shown bv lBeauge et al.l ( l2008h that orbital 
fits of multi-planetary systems in a mean-motion resonance (hereafter MMR) may be highly 
unreliable, though formal uncertainties of estimations may be apparently small. In these 
cases, the shape of the likelihood function may be complicated and may possess multiple 
comparable local maxima. Often this shape is model-dependent: addition of extra model 
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components leads to qualitative changes of the set of likelihood maxima. Often, every such 
local maximum provides a good fit of the RV curve but nevertheless is far from the actual 
orbital configuration. Such situation indicates one of the following items: 

1. The adopted model is imperfect. Some extra terms were not included or the terms in- 
cluded are wrong. 

2. The data are imperfect. The errors may have a non-Gaussian distribution, they may be 
correlated or they may incorporate some extra time-variable systematic part. Also, the 
data may cover too small time base or simply the number of observations is too small. 

Formally, we could use one of the local maxima (e.g. the global maximum) to construct an 
estimation of the parameters of the system. However, such estimation would appear strongly 
biased (usually to higher eccentricities) and its formal uncertainties would strongly under- 
estimate real errors of parameters. Therefore, in addition to the formal goodness-of-fit mea- 
sur e, we need some indi cator of statistical reliability of our fits. 

iBeauge et al.l ll2008h performed (for the system of HD82943) several fits with truncated 
RV datasets to explore the sensitivity of current orbital fits to future RV measurements. 
Unfortunately, this approach require too time-consuming computations. Here we need some 
simple and rapid (though perhaps quite rough) test of the 'statisti cal hea l th' of our orbital 
fits. For this goal, we use below the following approach. Recall (,Lehmanl . ll983l . § 6.4) that 
the asymptotic (A' oo) approximation to the variance-covariance matrix of the estimations 
of p is calculated as the inverse of the Fisher information matrix Q having elements 



J 



j=U=l "full.)/ "P" 



(8) 



When we deal with a well-conditioned situation, the likelihood function can be quadratically 
approximated in the vicinity of its maximum using the quadratic term °<: 5p^Q5p only. This 
is the asymptotic behaviour of the likelihood function for large time series (i.e.. A' 
The mentioned quadratic term approximates the multidimensional graph of the likelihood 
function by a paraboloidal hypersurface. Other terms in the expansion of the likelihood func- 
tion are insignificant and do not distort this shape essentially. However, when the number of 
observations is not sufficient for reliable determination of the parameters of the model, the 
matrix Q become ill-conditioned, and the extra terms easily produce distortions leading to 
multiple local maxima of the likelihood function. 

To assess the reliability of a particular fit we could calculate the condition number of 
the information matrix Q (i.e. the ratio of the biggest eigenvalue and the smallest one). 
The larger is this number, the lower is the reliability of the corresponding orbital fit. This 
condition number should be compared with the number of observations, because high-order 
terms in the expansion of the likelihood function tend to zero when grows. To demonstrate 
this, let us write down the full expression of the Hessian matrix of the usual function 
(speaking more accurately, of the weighted average of squared residuals): 

In this expression, the first term determines the asymptotic behaviour of the likelihood func- 
tion in the vicinity of its maximum. The second term (summation) reflects the non-linearity 
of the RV model and contains the residuals (v — /i), which average decreases as ff{\/\/N) 
when A' grows. This implies that the perturbing term in ^ grows according to the law 
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&'[\/N), whereas the Fisher matrix grows according to ff{N). Therefore, the relative mag- 
nitude of the second term in (|9]l is about 1 /\/N. Let us now transform Q to its diagonal form 
with eigenvalues E\ . . .Ed on the diagonal (we assume that E\ is the biggest eigenvalue and 
Ed is the smallest one). We may expect that the typical magnitude of the elements of the 
Fisher matrix is of the order of \/E\Ed and, therefore, the typical magnitude of the second 
term in ([9]l is about y/EiEd / \/N. The full matrix H is well approximated by its asymptotic 
representation — 2Q, if relative deviations of the corresponding eigenvalues are small. If this 
is so, the topology of the graph of the likelihood function should not be very sensitive to 
the perturbing terms. We can expect that absolute deviations of the eigenvalues should have 
similar magnitude about \jE\Ed/N. But the relative deviation of the smallest eigenvalue Ed 
is as large as \/E\ /Ed/ VN. Therefore, the condition number E\/Ed of the Fisher informa- 
tion matrix should not exceed the number of observations for a statistically robust orbital 
fit. This is simply a specification of the general fact, known from the numerical analysis, 
that the condition number is a measure of sensitivity of a matrix to small perturbations. 

There is a small clause. The elements of Q are measured in different physical units and 
it would be incorrect to use the condition number of Q itself (since it would depend on the 
choice of the measurement units). Instead, we may use the condition number of the scaled 
information matrix Q having elements Qij = Qij / ^/QuQjj- 

Unfortunately, the indicator depends on the epoch for which we obtain the fit, espe- 
cially when the RV model incorporates long-term trends. To obtain an informative estima- 
tion of we should choose the reference epoch near the middle of the time series span 
(or, when we merge datasets with different characteristics, near the weighted average of 
the timings r,j). In this position, the value of will be close to its minimum and, simul- 
taneously, the set of parameters will possess the best statistical properties (e.g., the mutual 
correlations of different parameters will be minimized). Even if absolute value of the quan- 
tity represents a too rough indicator of the fit robustness, this indicator seems to work 
quite well for the goal of intercomparison of different orbital fits for the same planetary 
system. To illustrate the use of the indicator we give here its (minimum) values for sev- 
eral planetary systems with well-determined orbits: "rf = 10 with A' = 409 (51Peg), "if = 6.7 
with N = 109 (70Vir), = 9.1 with N = 203 (14Her). Less perfect cases: ^ = 17 with 
N = 75 (HD69830), <^ = 120 with N = 487 (55Cnc). For the system of HD82943, we ob- 
tain ^ = 200 with A' = 165 (three-planet model with the eccentricity of the outermost planet 
fixed at zero). 



4 Best-fitting orbital solutions for HD37124 

From now on, we consider three main models of the RV data for HD37124. All models 
incorporate three common Keplerian terms as in the eq. the constant velocity terms 
(separate for the three datasets), and the annual term for the ELODIE dataset. The first RV 
model (I) does not incorporate anything else and thus has cl = 20 degrees of freedom. The 
second one (II) incorporates also an annual term for the Keck dataset and has d = 22. The 
third one (III) incorporates the same terms as (I) plus a linear trend (common for all the three 
datasets) and thus has d = 2\. For all these models, the ratio d/N kO.16 means that we are 
left with only w 6 observations per one parameter to be estimated. This indicates that the 
problem of obtaining a suitable orbital configuration of the system cannot be solved easily. 

Let us try to reconstruct the topological structure of the multidimensional graph of the 
likelihood function Since the dimension of the problem is large, we cannot look on 
corresponding hypersurface directly. Let us pick two (of d) free parameters x,y (say the 
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Fig. 3 Contour maps of the likelihood goodness-of-fit statistic min/ for RV fits of HD37124. The maps are 
plotted in the plane of orbital periods Pc and (measured in days). For each point in these panels, the 
likelihood goodness-of-fit function / was minimized over the rest of free parameters. Thick straight line mark 
the position of the 2/1 MMR. The level values of the plotted function are shown below the respective panels. 
The panels in different columns correspond to the different RV models (model 1, model 11, model III, from 
left to right). The top raw of panels corresponds to the fits based on the whole accessible data array, whereas 
the bottom raw corresponds to the fits based on the Keck data only. Note the differences in period ranges 
shown in top and bottom graphs. 



orbital periods Pc and P^i of the outer planets). Then we consider the function l'{x,y) = 
min/, where the minimization of the goodness-of-fit function I is performed over the rest 
of free parameters. This means that for any manually assigned values of x,y we obtain 
corresponding best-fitting values of other parameters and find corresponding goodness-of- 
fit measure /. Further, the resulting function of two variables can be visualised on a two- 
dimensional grid. 

Fig. [3] shows such plot in the plane of orbital periods Pc and P4. We can see that the 
likelihood function constructed from the full available dataset (top panels of Fig.[3]l has two 
main maxima with comparable values of I. The first one is centred on Pc « 2000 days and 
Pel w 870 days, and the second one on Pc w 1800 days and P^ w 900 days. We can see that the 
latter solution is close to the 2/1 MMR of the outer planets. Such orbital configurations are 
remarkable because only low-order MMRs can prevent planets on high-eccentricity orbits 
from close approaches and hence can make the whole system stable. This 'double-headed' 
shape of the likelihood function looks rather stable with respect to the choice of the RV 
model. From now on, we consider mainly the two mentioned families of solutions: the first 
one corresponds to the resonance 2/1 between planets 'c' and 'd', the second one is outside 
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Table 1 Best-fitting orbital solutions for the planetary system around HD 37124. 



parameter 


lA, 2/1 


IB 


IIA, 2/1 


IIB 


IIIA, 2/1 


IIIB 








planet b 








P [days] 


154.37(13) 


154.48(12) 


154.53(12) 


154.58(12) 


154.37(13) 


154.49(12) 


K [m/s] 


28.39(82) 


28.49(85) 


27.42(73) 


28.39(85) 


28.06(83) 


27.98(81) 


AH 


117.7(1.9) 


119.1(1.8) 


118.2(1.8) 


118.8(1.7) 


118.0(1.8) 


120.2(1.7) 


e 


0.116(33) 


0.116(33) 


0.108(29) 


0.081(27) 


0.093(33) 


0.064(28) 


to["] 


131(15) 


132(23) 


134(14) 


132(23) 


135(18) 


140(25) 


msin; [Mj„,,] 


0.635(46) 


0.637(47) 


0.613(44) 


0.635(46) 


0.628(46) 


0.626(46) 


a [AU] 


0.518(17) 


0.519(17) 


0.519(17) 


0.518(17) 


0.518(17) 


0.518(17) 








planet d 








P [days] 


912.5(7.7) 


867.7(9.1) 


902.3(8.7) 


869.3(9.5) 


904(14) 


868.3(8.5) 


^[m/s] 


13.8(1.2) 


13.3(1.0) 


15.1(1.0) 


13.70(97) 


13.1(1.1) 


14.1(1.0) 


A[°] 


290.3(7.4) 


345.9(4.3) 


289.5(5.4) 


344.7(4.2) 


297.8(8.3) 


346.0(4.0) 




U.jUjI^oo ) 


U.UjOl^ 12.) 


n A A 1 ( f^^ \ 
U.44HOJ ) 


n om/" Q'7^ 
U.ZU '[o') 


n Id/ ^ n\ 
O.jyylO) 


U.UZJ^Do ) 




283(12) 


i5{n) 


276(13) 


49(19) 


288(17) 


390(150) 


msin/ [Mj„/j] 


0.559(61) 


0.528(54) 


0.609(60) 


0.545(53) 


0.527(58) 


0.560(55) 


[AU] 






1 (^S9('S7'^ 

i .UoZl^ J / j 


1 .Ot i I JO ) 


i .UO'+^ J7 j 


I .Ojy ^ JO j 








plflJlCt c 










1782(41) 




1777(24) 


198nf55'l 


1776(41) 




K [m/s] 


14 9(] \) 


17.1(2.5) 


15.2(1.0) 


17.5(2.4) 




16.7(2.4) 


/L L J 


319 8f6 1 ) 


309 3(5 21 


322.5(4.5) 


312.7(4.6) 


311.9(7.1) 


306 2(5 71 


6 


yj.DZ. / y 




U. JDj JO ) 


U. J JJ \' ^) 


yj.H^O 1 yoj j 


U. J / J l^O'H- j 


W I \ 


LJ\j.z.\y .yj j 


325.8(7.6) 


251.9(6.6) 






QOn Qf'S 7'1 




I'^Adl) 

\J . / .^*T \ 1 1 } 


0.895(61) 


0.766(72) 


Ql 5^62^1 


720('76'l 




a [AU] 


2 648l'97^ 


2.85(11) 


2.643(91) 


2 941' 1 1"! 


2.642(97) 


2.80(11) 


^ '-s-yr-' 










0.58(33) 


0.78(27) 








ELODIE dataset 








Co [m/s] 


65.5(4.3) 


62.7(4.2) 


66.5(4.5) 


61.2(4.5) 


67.4(4.4) 


65.6(4.3) 


A [m/s] 


17.6(3.2) 


17.8(2.6) 


17.5(3.5) 


16.9(2.6) 


18.2(3.5) 


17.2(2.9) 


T [days] 


167(20) 


178(20) 


164(20) 


-181(22) 


162(18) 


171(20) 


fTjr [m/s] 


6.6(2.3 j 


5.9(2.3) 


7.4(2.3) 


6.5(2.3 j 


6.3(2.3 j 


j.6(2.4j 


r.m.s. [m/s] 


11.85 


10.96 


12.23 


11.21 


11.79 


10.45 






CORALIE dataset 








Co [m/s] 


4.3(4.2) 


7.9(4.2) 


4.8(4.3) 


7.4(4.2) 


5.3(4.4) 


8.6(4.4) 


(J* [m/s] 


11.3(4.0) 


10.9(4.0) 


12.3(4.2) 


11.2(4.1) 


12.1(4.1) 


12.0(4.1) 


r.m.s. [m/s] 


17.45 


16.93 


17.62 


16.73 


17.68 


17.38 








Keck dataset 








Co [m/s] 


7.56(90) 


9.11(84) 


6.8(1.0) 


10.8(1.1) 


7.47(80) 


9.00(79) 


A [m/s] 






4.5(1.4) 


3.8(1.6) 






T [days] 






-148(19) 


102(14) 






fJ* [m/s] 


2.58(54) 


2.48(55) 


1.24(71) 


2.02(56) 


2.44(54) 


2.01(57) 


r.m.s. [m/s] 


3.704 


3.713 


3.044 


3.381 


3.566 


3.359 






20 


22 


21 




f [m/s] 


8.269 


8.025 


7.816 


7.888 


8.180 


7.711 




213 


54 


190 


74 


145 


81 



The values of cq for ELODIE and CORALIE are given relatively to their first measuments. The uncertainties 
of the estimations are given in the brackets (e.g., 0.30(10) means 0.30±0.10, and 30.0(1.0) means 30. 0± 
1.0). The values for the mean longitudes X and time shift parameters T are given for the epoch 7D2452000. 
The uncertainties of the minimum masses msin; and of the semi-major axes a incorporate the 10% uncertainty 
of the stellar mass. The estimations of the effective RV jitters incorporate an analytic con'ection of the 
statistical bias as discussed in ( .Baluev... 2008c.) . 
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of this resonance (but may cover some other MMRs of relatively low, e.g. 7/3 and 5/2). 
Hereafter, we use the notation 'A' for the first family and 'B' for the second one. 

The full sets of estimated parameters for the three RV models are shown in Table [T] 
The respective minimum values of I depend on the model adopted. For the models I and III, 
the 'B' solution provides formally better fit to the RV data in comparison with the 'A' one, 
but for the model II the corresponding values of / are similar. Either model II or model III 
provide some improvement in the goodness-of-fit /, with respect to the model I. 

Let us note that the similar structure of the likelihood function can be also seen in the 
graphs constructed in the similar way for the RV model I, but with the use of the Keck data 
only (left-bottom panel in Fig.[3]l. However, the shape of the likelihood surface is more com- 
plicated in this case: the broad (much broader than in the top pictures) peak corresponding 
to the B-family is actually 'double-headed' itself (i.e., it consists of two close peaks having 
Pel w 850 days and P^i w 820 days). Moreover, the shape of the likelihood function con- 
structed with the use of the Keck data only is severely dependent on the adopted RV model. 
For the model II, the B-family peak is clearly split into two distinct and very distorted peaks 
(one with Pc « 1900 days, Pj w 870, and another one with Pc w 2200 days and P^ less than 
800 days, see middle-bottom panel in Fig. [3]l. For the model III, all the former peak are 
merged into a single very broad peak (formally centred not far from the the 2/1 MMR), 
which is continued to infinite P^ (right-bottom panel in Fig.[3j. Such behaviour, the further 
splitting of the local likelihood maxima and their severe sensitivity to the RV model, indi- 
cates that the analysis of the Keck data alone would yield sigruficantly less reliable results 
than the joint analysis of all available data. 

However, even with the use of the full RV dataset, almost all of the best fits possess 
large values of the condition number "rf, especially for the 'A' solution. This means that we 
should treat our results with a more care. In fact, no one of the best (in the sense of the 
goodness-of-fit measure) fits from Table [T] can be accepted. The values of the eccentricity 
Cc (and often those of the as well) are large and lead to a very soon disintegration of the 
corresponding orbital configurations. For comparison. Table [2] contains the estimations of 
parameters for the system with the eccentricity Cc or both the eccentricities 6^6^ fixed at 
zero. These fits have much smaller values of though worse goodness-of-fit. Numerical 
integration showed their dynamical stability and regular evolution on the time scales of (at 
least) 10^ yr, except for the fit II' A which showed some signs of chaoticity at the time scale 
of ~ 10^ yr, evidently due to a large = 0.337. It is important that only the 'A' group of 
solutions contains best-fitting circular orbits: the best-fitting solutions from the 'B' group 
approach the 2/1 MMR and softly turn into 'A group of solutions when 6^ decreases. 

From the interplay of the indicators and / described above, we derive the following 
conclusion. Although the orbital fits from Table[T]show relatively small scatter of the resid- 
uals, this small scattering is in fact fictitious, as indicated by the corresponding values of "rf. 
The number of RV data points and their temporal coverage still cannot constrain the full set 
in the model parameters reliably. The full RV model is 'overloaded'. Injecting some kind 
of a priori information (e.g., fixing the eccentricities at low values) allows to overcome the 
obstacle of the statistical ill-determinacy. The resulting orbital fits possess better statistical 
reliability, but by the cost of some increase of the RV residuals scatter. Nevertheless, this 
increase is necessary to obtain a physically realistic orbital configuration. 

However, we cannot rule out the possibility that the actual orbits of the planets 'c','d' 
are far from circular. To find more realistic orbital configurations than those from Table [T] 
but corresponding to eccentric orbits, we need to account for more subtle requirements of 
the dynamical stability in our analysis. 
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Table 2 Low-eccentricity orbital solutions for the planetary system around HD 37124. 

parameter I' A, 2/1 I" A, 2/1 ll'A, 2/1 II" A, 2/1 III'A, 2/1 III" A, 2/1 



planet b 



r [uaysj 


1 S4 1 9"! 
1 J^l^ IZJ 


1 S4 'KA ( M\ 


1 Jt. JZl^ IZJ 


1 S4 1 '^'1 
1 J't. JOl^ ID j 


1 S4 ^^M 1 "1 
1 J4. J J 1. 1 1 J 




K [m/s] 


n o o 1 

2o.5l(y 7) 




27.60(92) 


29.32(9 /j 


27.96(83) 


25.20(82) 


A L J 


\ \ n A (1 t\\ 

11/ .^l^Z.U J 


117 1 1 Ci\ 
li 1 .lyl.y) 


i i / .L)( i .oj 


1 1 7 rif" 1 Ci\ 

i i 1 X}\ i .VJ 


1 1 Q ^f"! 1\ 
I lo.Ol^ ^ ■ 1 j 


1 1 y ^ 1 7^ 

iio.J(i./J 


c 






u.uyu^zo J 


n HQ/i /"i 1 ^ 
U.Uo4(jl ) 


U.UOz^zy ) 


U. Do J (29 ) 


to L J 




111 /ti\ 
i J i (^zJ j 


1 1 A/ 1 Q^ 
1 lo^ivj 


124(22 j 


1 At^("l^\ 
140(^ZJ ) 


14/(24 J 


msin; [M/„,,] 


n e^A A i AQ'\ 
U. 044(^45 J 




U.Oi / (40 ) 


U.OJO(4fS J 


fl /^Tif AA\ 
U.OZj^4oJ 


U.oJl (4o ) 


a [AU] 


0.518(17)5 


0.518(17)4 


0.518(17)9 


0.518(17)5 


0.518(17)5 


0.518(17)5 








planet d 








P [days] 


908(11). 


905(11). 8 


899.3(6.9) 


904(1 2). 1 


883(10). 7 


881.9(8.5) 


K[m/s] 


15.2(1.1) 


14.17(92) 


17.1(1.2) 


14.13(94) 


15.65(88) 


15.21(81) 


xn 


322.3(4.1) 


324.3(4.4) 


322.9(3.2) 


324.3(4.3) 


328.9(3.6) 


329.8(3.6) 


e 


0.174(88) 


O(fixed) 


0.337(71) 


O(fixed) 


0.095(71) 


O(fixed) 


con 


349(21) 




1.7(9.8) 


- 


358(33) 


- 


msini [Mj,,,,] 


0.614(62) 


0.572(53) 


0.689(66) 


0.570(54) 


0.626(54) 


0.608(52) 


a [AU] 


1.690(58)1 


1.687(58)4 


1.679(57)3 


1.685(58)2 


1.659(57)8 


1.657(56)6 








planet c 








P [days] 


1810(41).3 


1839(50).0 


1759(28).0 


1834(54) .7 


1815(40).0 


1822(44) .7 


K [m/s] 


13.4(1.3) 


12.0(1.1) 


15.4(1.5) 


11.9(1.0) 


12.5(1.1) 


11.86(88) 


A [°] 


310.7(4.7) 


308.5(5.2) 


318.8(3.7) 


310.5(5.2) 


303.0(5.1) 


301.8(4.7) 


msin/ [M;„p] 


0.682(81) 


0.611(67) 


0.774(89) 


0.607(67) 


0.638(68) 


0.604(60) 


a [AU] 


2.678(98)0 


2.70(10)61 


2.627(92)3 


2.70(10)19 


2.682(98)5 


2.690(99)1 


'^l 


- 






- 


1.18(26) 


1.22(24) 






ELODIE dataset 








Co [m/s] 


t^A Q\ 
04.Z(4.0 ) 


64.9(4.3) 


63.9(4.5) 


DJ.U(^4.4j 


e^Q A ( A ^ \ 
05.4(^4. J ) 


by. i (^4.3 ) 


A [m/s] 


18.9(3.1) 


19.9(3.3) 


17.8(3.3) 


19.8(3.4) 


18.9(3.7) 


19.5(3.7) 


T [days] 


167(18) 


164(17) 


167(20) 


164(17) 


158(17) 


156(16) 


OV [m/s] 


6.8(2.2) 


7.1(2.2) 


7.6(2.2) 


7.4(2.2) 


6.7(2.3) 


6.7(2.2) 


r.m.s. [m/s] 


11.95 


12.33 


12.53 


12.45 


12.08 


12.18 






CORALIE dataset 








Co [m/s] 


5.0(4.2) 


4.5(4.0) 


5.1(4.5) 


4.2(4.0) 


6.5(4.5) 


6.3(4.3) 


(J* [m/s] 


11.6(4.0) 


10.5(3.8) 


12.8(4.2) 


10.3(3.9) 


12.7(4.1) 


12.1(4.0) 


r.m.s. [m/s] 


17.42 


16.76 


17.68 


16.53 


17.67 


17.35 








Keck dataset 








Co [m/s] 


7.56(75) 


7.89(76) 


9.52(93) 


8.9(1.0) 


7.88(61) 


8.04(61) 


A [m/s] 






4.9(1.3) 


2.0(1.4) 






T [days] 






134(13) 


103(34) 






OV [m/s] 


3.36(55) 


3.51(56) 


2.21(55) 


3.41(55) 


2.15(56) 


2.21(57) 


r.m.s. [m/s] 


4.218 


4.319 


3.436 


4.193 


3.394 


4.446 



d 18 16 20 18 19 17 

/"[m/s] 8.710 8.776 8.280 8.769 8.065 8.057 

34 29 40 32 36 31 

The same notes as in Table[lJto be applied here. In these fits, the eccentricity or both the eccentricities 
£(. and ej were fixed at zero. All these configurations con'espond to the 2/1 MMR. The estimations of the 
semi-major axes and orbital periods are given with one or two excessive decimal digits (shown after two- 
digit uncertainties enclosed in brackets) to allow an unambigious reproducing of the long-term dynamics of 
these configurations (see Section [xTt . For example, 1810(41). 3 means 1810. 3±41 and 2.70(10)19 means 
2.7019±0.10. 
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5 Dynamical interpretation 

To obtain more realistic stable orbital configurations for this planetary system, we continue 
to use the method of planar plots of partially minimized goodness-of-fit statistic I. But now 
we examine orbital solution from a two-dimensional grid more carefully: for each solution, 
we perform a numerical integration in order to rule out rapidly disintegrating configurations. 
To perform such integrations, we need to know true masses of planets in the system. As it 
can be seen fror n llHi, they depend on the mass of the star and on the orbital inclina- 
tions. Following IVogt et all l l2005b . we adopt = O.TSMs with an uncertainty of 10%. 
Unfortunately, using the Keplerian RV model we can estimate only the minimum masses 
m sin i where the inclination i remains unknown. Until the gravitational interactions between 
planets in the system are directly observed in the RV curve, the best thing that we can do is 
to assume a priori that the orbits are coplanar with / = 90° . If the actual inclination is less 
than 90°, the true masses of the interacting planets become larger and the stability region of 
the system become more narrow than for the edge-on configurations. The same effect is ex- 
pected from non-zero mutual inclinations of the orbits, due to the well-known phenomenon 
of the e-i coupling. 

Since much troubles in obtaining a realistic orbital configuration of the system are due 
to the large eccentricity of the outermost planet, let us firstly consider the plane of eccentric 
variables {ecCosC0c,ecS'mC0c}. The corresponding maps are plotted in Fig.|4l We can see 
clearly the sophisticated shape of the likelihood surface: among the 'A' families of solu- 
tions, no one possess a single maximum. Instead, we can see 2-3 local maxima; all of them 
correspond to large values of Cc- The 'B' family shows single maximum, but again at large 
Cc- No one of these local maxima corresponds to a stable configuration. Stable solutions oc- 
cupy only regions of small or moderate Cc- It is important to note that the small-eccentricity 
solutions correspond to the A' configuration only: when Cc decreases, a B-type solution 
approach the 2/1 resonance and softly turns into an A-type one. 

As it can be seen from Fig. |4l we can easily find stable solutions from the A' family. 
Such solutions can possess values of / as small as w 8.6 m/s (model I), w 8.2 m/s (model II), 
and w 8.1 m/s (model III). It is interesting that the region of stable configurations is some- 
what correlated with one of the local minima of I in the 'A' layer, located in the third quadrant 
of the coordinate plane. On contrary, we face much difficulties with obtaining a stable con- 
figuration from the family 'B'. The main reason for almost all 'B' solutions to be unstable 
is that the best-fitting period ratio Pc/Pd « 2.1 — 2.3 is quite small and is not fixed (with a 
necessary precision) at any MMR of low order. 

To locate stable solutions from the 'B' group, we use another pair of variables. In Fig.|5] 
the partially minimized I is plotted in the plane {Pc,ec). We can see that the fits with low 
Cc and with Pc fixed far from the 2/1 resonance, possess uncomfortably large values of 
/. Fig [6] illustrates this more clearly. In this figure, we plot the function /, minimized so 
that the eccentricity was fixed at a safe value of 0.15 and the period ratio Pc/Pd was 
fixed at the values marked on the abscissas. Since the value of the eccentricity was fixed 
at a small value, the effects of the RV model non-linearity are significantly suppressed (it 
follows from relatively small values of in Table O, and thus we may try to find some 
confidence intervals for the ratio Pc/Pd using the standard asymptotic (A' — > °°) theory of 
point estimations. We can see that the values Pc/Pd > 2.5 and Pc/Pj < 1.87 are beyond 
the 99% confidence interval, when we use the RV model I. For the other RV models, this 
confidence interval becomes even smaller: Pc/Pd G [1 .86, 2. 1 1] for the model II and Pc/Pd G 
[1.95,2.3] for the model III. When the eccentricity decreases, these confidence regions 
are shrinking around Pc/Pd = 2. 



14 



0.6 

0.4 

a' 0.2 
o> 
<u 
E 

I ° 
o 

J -0.2 
-0.4 
-0.6 



9.17 

9.00 

8.83 

8.67 

8.51 

3.35 

3.27 



0.6 

0.4 

3 0.2 
E 

a)'-0.2 
-0.4 
-0.6 



, / 



9.07 • 
8.90 • 
8.73 ■ 
3.57 ■ 
8.41 ■ 
8.26 ■ 
8.11 
8.03 ■ 



-0.6 -0.4 -0.2 0.2 0.4 0.6 



-0.4 -0.2 0.2 0.4 0.6 



0.6 
0.4 



E 

& 



m' -0.2 
-0.4 
-0.6 



1 1 ! 1 1 1 1 










N 






















' i 









































8.85 - - - 

8.68 - ■ ■ 

8.52 

8.36 

8.20 

8.05 

7.89 

7.82 



-0.6 -0.4 -0.2 



0.2 0.4 0.6 



co' 0.2 

g" 

E 

t ° 
«' -0.2 

-0.4 

-0.6 





















\ 


















' .' / 










' i 










■ \ v. \ 
'■ \ \ 










1 









9.00 ■ 

8.67 ■ 

8.50 ■ 

8.34 ■ 

8.19 • 

8.04 ■ 

7.96 • 



-0.4 
-0.6 



0.6 
0.4 













1 1 


- 


♦\ 


v...'' 

) 




















































, y 1 

; i 


D 






















-0 


6 -0.4 -0.2 


3 0.2 0.4 0.6 
















•-".y..- 






























\ '«\ 












'\\\\ 










: ; 










1 


\ ■■ ' 

\_ 

















8.93 ■ 

8.76 • 

8.60 ■ 

8.43 ■ 

8.27 ■ 

8.12 ■ 
7.97 



8.73 ■ 

8.56 ■ 

8.40 ■ 

8.25 ■ 

8.09 ■ 

7.94 • 
7.79 

7.72 ■ 



-0.6 -0.4 -0.2 0.2 0.4 0.6 
e_c cos(omega_c) 



-0.6 -0.4 -0.2 0.2 0.4 0.6 
e_c cos(omega_c) 



Fig. 4 Contour maps of the likelihood goodness-of-fit statistic for RV fits of HD37124. These maps are 
plotted on the plane {e^ cos (Oc,ec sinWc) in the same way as in Fig. [3] Solutions which coiTespond to orbital 
configurations disintegrating in less than 30000 years are mai'ked by bold grey dots, other solutions are 
marked by fine (green in the electronic version of the paper) dots. These dots are arranged according to the 
polar coordinate system. Plots in the left column correspond to the 'A' solutions close to the 2/1 resonance 
(with no more than 5% relative deviation of the period ratio), plots in the right column correspond to the 'B' 
solutions. The white regions (where < 0.7) mark the points for which the fitting algorithm could not find a 
solution from the coiTesponding family and switched to another one (that appeared significantly more likely). 
The 'A' and 'B' families of solutions overlap in the region of large without possibility of any smooth seam, 
but they seem to be sewed smoothly in the region of small e^- The top, middle, and bottom pairs of panels 
correspond to the RV models 1, II, and III, respectively. 
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Fig. 5 Contour maps of 
the likelihood goodness-of- 
fit statistic for RV fits of 
HD37124. These maps are 
plotted on the plane {P^ec) 
in the same way as in Figs. [3] 
and \4\ Three panels corre- 
spond to the models 1, II, 
and III (from top to bot- 
tom). In each of these pan- 
els, the A and B families of 
solutions were merged in a 
single plot. The bold nearly 
vertical lines mark the so- 
lutions having ratio of the 
best fitting periods and 
Pj close to the 2/1, 7/3, 5/2, 
and 8/3 commensurabilities 
(lines from left to light). 
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Fig. 6 Graphs of the likelihood goodness-of-fit function /, which was minimized so that the eccentricity 
was fixed at 0.15 and the ratio of orbital periods Pc/Pj was fixed at the values shown on the abscissas. Three 
panels correspond to the RV model I, II, and III (from left to right). The bold horizontal lines shows the levels 
of min/ yielding the 80%, 95%, and 99% asymptotic confidence intervals for Pc/Pd- These levels correspond 
to the values of the likeli hood ratio statis tic, which provide the asymptotic false alarm probabilities as small 
as 20%, 5%, and 1% (see lBaluevl , r2008d) . 
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However, we can note a promising region of high-eccentricity solutions near the reso- 
nance 5/2 of the two outer planets. Surprisingly, numerical integration of the best fitting 
configuration with Pc fixed at 2170 days and fixed at 0.4 (in the RV model I) showed 
a quite regular evolution without any signs of instability at the time scale of at least 10^ 
years. The evolution of this orbital configuration is a large-amplitude oscillation around an 
antialigned apsidal corotation. The refore, this solution belong s to the class of orbital config- 
urations in the 5/2 MMR found in ( Gozd ziewski et al. l|2006^ using only the Keck data. 

It seems that, due to the non-linearity of the RV models coupled with lack of the data 
and insufficient time coverage, none of the local minima of / lies near the real orbital config- 
uration of this system. Probably, these multiple local minima are only the fictitious 'ripples' 
produced by the lack of the observations. All these 'ripples' are located deeply in the zone of 
dynamical instability, and therefore cannot be close to the real configuration of the system. 
It is possible to find strictly the best-fitting orbital solution(s) simultaneously satisfying the 
stability requirement. Evidently, such solutions would be attracted by one of the 'ripples', 
and therefore would be close to the boundary of the domain of system stability. Therefore, 
we would have to deal with large difficulties concerni ng the very complicated s tructure 



of the parameter space near such boundarie s (see, e.g., Gozdziewski et al.l l2008h . which 
probably represents the Arnold web (see, e.g. lFroeschl6 et all 200d) . Near such boundaries, 
the dynamics of the system is very sensitive to small changes of parameters, the stability 
map of the parametric space is strongly irregular, and hence the resolution of the parame- 
ter space should be chosen dense enough. This requires very time-consuming calculations 
for checking the stability of probe orbital configurations in this region. On contrary, it does 
not look likely that real planetary systems can be found in such extremely dynamically ac- 
tive regions: it would be rather difficult to explain how the system could migrate (without 
disintegration) to such a state through the dense web of the instability threads and why it 
stopped in a thin island of stability instead of moving further to dynamically unstable con- 
figurations. Therefore, the reliability of such 'hardly-stable' solutions would be too low to 
justify the associated time-consuming calculations. In addition, considering the boundary of 
the stability domain, it is rather difficult to understand the physical mechanism stabilizing 
the configurations in the gi ven domain. 

In this paper, following iHadiidemetrioul bOOd) , we will pay more attention to the cen- 
tres of the stability domains which point out families of orbital configurations having some 
'stock of stability'. For having a clear picture of possible dynamical regimes of a planetary 
system, it is the position of the centre of the stability domain which should be located and 
for which we should know possible uncertainties. Considering the kernel of a stability do- 
main, we avoid dealing with the sophisticated dynamical structure near its boundary, which 
is hardly able to carry much information about the dynamics of the real system. When an 
orbital solution has some stock of stability, the dynamics of the corresponding planetary 
configuration is much more regular and much less sensitive to small changes of parameters. 
To obtain such orbital solutions, we will try to decrease the dimension of the problem (i.e., 
the number of degrees of freedom d) using certain a priori information about the stability 
of resonance planetary systems. 



6 The value of apsidal corotation resonances 

Bearin g in mind the results bv lji et al.l ( l2003h : lHadiidemetrioj ( l200^ . l2008h : l\^atzis & Hadiidemetrioul 
( l2006h . let us recall that regular stable motions on high-eccentricity orbits with small period 
ratio are only possible if the planets are trapped in a MMR and simultaneously are close 
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to (or, at least, not far from) an apsidal corotaion resonance. The details of the theory of 
apsidal corotation reso nances, along with necessary formulae and further references can be 
found, for instance, in teeauge et all boOS). For a brief summary, let us consider two planets 
trapped in the p/q MMR, i.e. having the ratio of orbital periods ~ p/q with p > q. 

We can write down the resonant angles 

ph - qh ph-qh , . „, 

Sl = COi, S2 = 0)2 (10) 

p-q p-q 

and the canonically conjugated action variables 



/i =Li (^l-^l-e2j, i^=L2(l-^Jl-elj, (11) 

where A,- are the mean longitudes of the planets and L,- :^ m, sfa^ are the Delaunay action 
variables. After averaging the Hamiltonian H of the system over the fast variables (i.e, over 
the mean longitudes A,,) keeping the resting slow ones, the resulting averaged Hamiltonian 
(//) depends on the canonical variables i,-,/,, and (as on parameters) on the masses m,- of the 
planets. Evidently, this averaging accounts properly for the orbital resonance. The averaged 
equations of motion are then given by 

dt dsi ' dt dli 

Suppose that some values s* , I* determine the position of an extremum of (H) . We can easily 
see that every such extremum provides a stationary solution = s* , /,■ = /* of the averaged 
system i ll2b . Such stationary solution is often called 'apsidal corotation resonance' (here- 
after ACR). If the initial state of the planetary system slightly deviates from an exact ACR, 
the motion is a stable oscillation around the exact stationary solution, because the planets 
are prevented from close approaches. If the orbits of planets are highly eccentric and are far 
from stationary solutions of the averaged Hamiltonian equations, the motion is, most proba- 
bly, highly chaotic and unstable: the secular drift of resonant angles jlOl l leads the planets to 
too close approaches destabilizing the system. Stable solutions with one or both the resonant 
angles circulating are also possible in some cases; however, for high-eccentricity configura- 
tions, the ACRs mark centres of dynamical stability (see e.g. Hadjidem etriou . 2006 , 2 00^) . 

To obtain nominal orbital configurations of the system, we require from the resonant 
planets 'c' and 'd' to be locked in an exact ACR (while neglecting the influenc e of the irmer- 
most p lanet 'b'). This can be justified not only by the stability considerations. iBeauge et al] 
showed that adiabatic dissipative perturbations (e.g., interaction with a protoplane- 
tary disk) can cause planet pairs in a MMR to be captured in an ACR as well. 

The requirement of the ACR lock implies four algebraic equations: d{H) /dscj = and 
d{H) /died = 0. We neglect here the gravitational influence of the innermost planet 'b': it is 
seemingly non-resonant with the outer planets and probably should not affect their resonant 
dynamics much. The four equations mentioned above put certain constraints on the full set of 
free parameters to be estimated from RV data and decrease the number of degrees of freedom 
by four. It is very important, because this decreasing makes the problem significantly better 
determined: we have about 8 observations per a degree of freedom instead of about 6. 

More detailed description of the algorithm used to obtain such ACR fits, is given in Ap- 
pendix |A] The resulting fits of ACR solutions are given in Table [3] It is important that since 
the position of the ACR depends on the planetary mass ratio only and is almost independent 
of the individual planet masses, the solutions from Table [3] are almost independent of the 
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Table 3 ACR fits for the planetary system around HD 37124. 



parameter 


FA, 2/1 


FB, 5/2 


IF A, 2/1 


IFB, 5/2 


UFA, 2/1 


IIFB, 5/2 








planet b 








r [daysj 




lD4.jo(12j 


i j4.40( 12) 


1 j4. jo(^ 13 ) 


1D4.38(^ 12) 


1^/1 izf ^ '^\ 
154.35(12) 


k [m/s] 


29.21(92) 


28.9(1.0) 


28.84(83) 


28.9(1.0) 


28.69(85) 


28.8(1.0) 


A [°] 


1 1 r\ o / 1 fi\ 

iiy.j(i.y) 


1 19.1(1.8) 


1 18.8(1.7) 


11A ^ ( ^ c\\ 
1 19.1(1.9) 


1 1 A 1 /" 1 O \ 

120.1(1.8) 


119.8(1.8) 


e 


().u/y(iz) 


r\ m /I / T /I \ 

U.U/4(34) 


0.095(29) 


A AT A / 'J A\ 

O.U/4(34) 


0.066(31) 


0.067(34) 


to ["] 


1 jU(zj) 


154(23) 


1 19(18) 


154(24) 


136(26) 


164(26) 


msin; [M/„,,] 


0.653 (4s) 


f\ £1 A'^ 1 Aa \ 

U.o47(4o) 


0.645(4 7) 


A £i A1 f Ar\\ 

0.647(49) 


A £1 A'^ f An\ 

0.642(4/) 


{\ £i A A 1 A \ 

0.644(48) 


a [AU] 


0.518( 17)6 


0.518(17)5 


0.518( 17)8 


0.518(17)5 


0.518( 17)6 


0.518(17)5 








planet d 








r [days J 


yu/.o(o.y J 


o /4.0^ /.8) 


8y /.2( /.8) 


8 /4. 4(^8. 2 ) 


896(1 1 ).i 


on A t^(n z\ 
8 /4.6( /.5 ) 


K [m/s] 


16.10(81) 


13.82(86) 


16.53(81) 


13.84(91) 


16.63(78) 


14.40(91) 


A 


TIC n n\ 

ji5.y(j.o) 


340.7(3.6) 


316.4(2.5 ) 


340.7(3. /) 


31 /. 4(2.9) 


339.7(3.5 ) 


e 


0.306(88) 


0.221(26) 


0.311(77) 


0.221(26) 


0.289(73) 


0.199(33) 


co["] 


320.8(8.8) 


154.5(8.4) 


334.4(7.9) 


154.6(8.5) 


315(11) 


154.0(8.6) 




0.650(54) 


551 (50"! 


665f55'l 


552f51 "1 


0.668(55) 


574(53'! 


a [AU] 


1.68y(57)7 


1.648(56)3 


1.676(57)8 


1.648(56)1 


1.675(58)6 


1.648(56)4 








planet c 








P [days] 


1815(17).3 


2186(20).4 


1794(15).5 


2186(21).0 


1792(22).8 


2187(19).0 


K[mlf.] 


12.6(1.1) 


10.54(93) 


12.8(1.1) 


10.54(95) 


12.5(1.0) 


10.30(87) 


X ["] 


311.8(5.3) 


301.0(5.5) 


318.2(5.0) 


301.0(5.6) 


307.7(6.0) 


300.3(5.5) 


c 


0.122(56) 


0.377(74) 


0.132(42) 


0.379(78) 


0.136(54) 


0.333(76) 


ffl ["] 


267(23) 


334.5(8.4) 


278(21) 


334.6(8.5) 


250(20) 


334.0(8.6) 


msin/ [Mj,,,,] 


0.63y(70) 


0.570(63) 


0.648(71) 


0.570(64) 


0.631(68) 


0.558(60) 


[AU] 


2.683(91)0 


3.03(10)69 


2.662(90)4 


3.03(10)66 


2.660(91)8 


3.03(10)75 




- 






- 


0.87(32) 


0.44(26) 






ELODIE dataset 








Co [m/s] 


64.7(4.4) 


65.5(4.7) 


64.6(4.6) 


65.5(4.5) 


67.4(4.5) 


66.8(4.5) 


A [m/s] 


18.8(3.1) 


20.2(3.2) 


18.1(3.3) 


'^A / T O \ 

20.2(3.3) 


1 O /I / T A \ 

18.4(3.4) 


1 A O / T O \ 

19.8(3.3) 


T [days] 


169(19) 


167(17) 


168(20) 


167(18) 


163(18) 


165(18) 


fJ* [m/s] 


7.2(2.2) 


7.7(2.1) 


8.2(2.1) 


7.7(2.1) 


6.8(2.2) 


7.6(2.1) 


r.m.s. [m/s] 


12.20 


12.60 


12.85 


12.51 


12.15 


12.33 






CORALIE dataset 








Co [m/s] 


4.8(4.1) 


6.2(3.8) 


4.5(4.1) 


6.1(3.8) 


5.7(4.4) 


6.3(3.9) 


/-T rm/cl 
\J-^ L 1 1 u ft J 


11 5(3 9) 


9.4(3.8) 


11.3(4.0) 


9.4(3.8) 


12 6(4 1) 


10 on 8"! 


r.m.s. [m/s] 


17.91 


15.91 


17.43 


15.82 


18.44 


16.06 








Keck dataset 








Co [m/s] 


7.13(73) 


8.13(72) 


9.1(1.0) 


8.1(1.0) 


7.11(66) 


8.15(71) 


A [m/s] 






3.9(1.3) 


0.1(1.1) 






T [days] 






131(17) 








ov [m/s] 


3.48(56) 


3.15(56) 


2.77(54) 


3.27(55) 


2.93(55) 


3.08(55) 


r.m.s. [m/s] 


4.322 


4.108 


3.762 


4.148 


3.910 


4.001 






16 


18 


17 




I [m/s] 


8.851 


8.616 


8.601 


8.693 


8.530 


8.549 




39 


56 


38 


64 


48 


79 



The same notes as in Table[lJto be applied here. These orbital elements were obtained for the Jacobi coordi- 
nate system with appropriate masses assigned to the reference barycentres. See Appendix|A]for more details 
concerning the algorithm used to obtain these fits. The semi-major axes and some orbital periods are given 
with excessive decimal digits, as in Table[2] 



19 



o 
< 



o 

CO 




-180 -120 -60 60 120 180 
RESONANT ANGLE s_d 



-180 -120 -60 60 120 180 
RESONANT ANGLE s_d 



Fig. 7 Top panels: orbits of the HD37124 system for tlie ACR fits (model 1 only). Left panel is for the 2/1 
MMR; right panel is for the 5/2 MMR. Straight solid lines mark the conjunction positions. Broken lines mark 
the lines of apses. Middle panels: contour maps of the averaged Hamiltonian {H^j) in the planes of resonant 
angles for the same MMRs and eccentricities as in top panels. Crosses mark the actual ACR positions 

con'esponding to local maxima of (Hcd)- 



assumptions about the orbital inclination of the system: (mesini)/(m^sini) = mc/nij. We 
only need to adjust the ratio of orbital periods by the quantity ^(to/M^,) (see formula ( 122b 
in Appendix [Ajl, but this adjustment have insignificant effect on the RV fit quality (until 
we consider inclinations as small as a few degrees). This invariance with respect to orbital 
inclination represents an extra advantage of the use of the ACR fitting procedure. 

We can see that in the case of the 2/1 resonance, the best fitting ACR is asymmetric 
with difference between the longitudes of periastra about 60°, whereas in the case of the 
5/2 resonance, the corresponding ACR is symmetric and antialigned (i.e., Oc — O);/ = 180°), 
see Fig. |7] The corresponding RV curves fit all available RV data satisfactorily, including 
the ELODIE measurements (Fig. The values of the goodness-of-fit measure / are not 
increased very much with respect to those from Table [T] and are comparable with those 
from Table 121 This means that one of the ACR fits can reflect the true configuration of this 
system quite well. The corresponding condition numbers for the 2/1 resonance in Table[3] 
are much less than in Table [T] This indicates that the topology of the likelihood surface 
becomes simpler and closer to the desirable paraboloidal one. The troubles connected with 
multiple local maxima of the likelihood have been overcome in the ACR fits. We can now 
say more definitely, that the effect of the putative annual term in the Keck RV data on the 
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Fig. 8 RV curves for the 
ACR solutions FA (top 
panel) and I'B (bottom 
panel), plotted together with 
the RV measurements and 
their residuals. The error 
bars do not incorporate 
estimated values of the 
RV jitter. The best-fitting 
annual sinusoidal drift of the 
ELODIE measurements was 
preliminarily subtracted. 
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fit quality is similar to the effect of the linear RV trend (which could be due to a long- 
period planet or brown dwarf in the system). These extra terms are significant in the fits 
for the 2/1 resonance (the 'A' group of solutions). For these fits, the estimations of false 
alarm p robabilities, calculated from the corresponding likelihood ratios according to lBaluevI 
( l2008ch , are about 2% and 0.3% (for the annual term and for the linear drift, respectively). 
The 5/2 resonance (B) solution does not require these terms in the RV model. This dilemma 
can be solved by future observations only. 



7 Long-term dynamics 

Now we consider the dynamical regimes of our nominal orbital configurations more closely. 
We are especially interested in how much the innermost planet can disturb the apsidal coro- 
tation of the two outer planets. This effect may be split in two categories: 

1. The planet 'b' can inspire some extra oscillation of the outer planets 'c', 'd' near their 
unperturbed ACR ('unperturbed' means obtained without taking into account the influ- 
ence of the planet 'b'). 

2. The planet 'b' can shift the position of the libration centre from the unperturbed ACR. 

Only the first effect may significantly affect the system stability, because only a large- 
amplitude oscillation around the libration centre can significantly increase the probability 
of close approaches of the resonant planets. 

Fig. m illustrates both effects. We can see that the orbital configuration FA taken from 
Table [3] shows moderate oscillation of the eccentricities ec,e^ and libration of resonant an- 
gles Sc,sj around some equilibrium values. However, the centres of oscillations are some- 
what shifted with respect to the unperturbed ACR. Heuristically, to counterbalance the grav- 
itational influence of the perturbing planet 'b', we need to increase somewhat the mass of 
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Fig. 9 Temporal evolution of the ACR configurations of planets in HD37124 (for RV model I only). Left 
panels siiow the evolution of the eccentricities, right panels show the evolution of resonant angles. Top pair: 
resonance 2/1, orbital parameters were taken from Table[3] solution I' A. Middle pair: the same case, but the 
value of the mass m^- was increased so that the value of increased by 1 m/s. Bottom pair: resonance 5/2 
(Tabled solution I' B). The same character of motion is conserved on the time scale of 10* yr (and probably 
much longer) for all the cases. 



the outermost planet 'c' . This assumption is confirmed by numerical integration: orbital con- 
figuration with iric increased by about 8% shows much less libration amplidute. Moreover, 
this adjustment decreases the statistic / by about 0.1 m/s (i.e., the scatter of the data around 
the RV model becomes slightly less). In the case of the 5/2 resonance, (solution FB), the 
libration amplitude is quite small for the unperturbed ACR solution already. The evolution 
of these orbital configurations appears perfectly regular and does not show any instability at 
the time scale of at least 10^ yr. 

More detailed analysis shows a great diversity of dynamical behaviour in the vicinity 
of the ACR solutions. The amplitude and character of the librations can be different. For 
certain orbital configurations inside the 2/1 MMR, the system may switch (from time to 
time) between alternating asymmetric ACRs with cOc — coj « —60° and cOc — (Od ~ +60°. 
Librations surrounding simultaneously the pair of stable asymmetric stationary solutions 
and unstable symmetric aligned ACR (see the left-middle panel in Fig.[7ll are also possible. 
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Fig. 10 The likelihood ratio periodograms of the RV residuals for different three-planet orbital models of the 
system of HD3 7 1 24. Every value of these periodograms represents the modified likelihood ratio statistic (see 
iBaluevl . [20080 °), calculated for the base model of the residuals (free constant velocity offsets + free common 
linear drift) and for the alternative model incorporating also a sinusoidal variation with free amplitude and 
phase. The panels in the left column represent the periodograms constructed from the residuals of full RV 
dataset. They were constructed in the range of periods starting from 10 days, since the errors in the dates 
of ELODIE measurements do not allow to fit accurately more short periods. The graph in the right columns 
show the periodograms constructed for the Keck residuals only. The top pair of panels is for the ACR solution 
in the 2/1 MIVIR (RV model I). The bottom pair of panels is for the ACR solution in the 5/2 MIVIR (also RV 
model I). For the case of graphs to the left, the normalized frequency bandwidth W ~ 350, and for the graphs 
to the right W « 3400. 



8 Testing the existence of extra planets 

When the best fitting orbital structure of the planetary system appears unstable, we can 
suspect that more planets orbit the star. It was the case for the system of /i Ara: the de- 
termination of rea l istic o rbits of jU Ara b and c represented an essential difficulty (e.g., 
Oozdziewski et al. ','2005') until the discovery of the planet ji Ara e placed all in the places 
JPepe et al.. .2007 ; Gozdziewski et al.L l2007h . To test the hypothesis of an extra planet or- 
biting HD37124 , we use the likelihood ratio-based periodogram describe d in the pape r 
iBaluev, '2008c'). This periodogram represent a generalization of the usual i Lombl l ll976h - 
[Scargle ( 1982) periodogram (as well as of the data-compensated discrete Fourier transform 
periodogram by iFerraz-Melld ( Il98lh ) and incorporates a built-in estimation of the RV jitter. 

The graphs of such periodograms of RV residuals for several orbital solutions and in- 
volving different datasets are shown in Fig. [TOl We can see that none of the peaks rises 
clearly beyond the apparent noise level. The formal significances of a periodogram peak can 
be assessed using the analytic expression for the associated false alarm probabilities (i.e., 
the probability to claim that the pe ak is sta t isticall y significant when actually it is a result of 
noise fluctuations) from the paper ( lBaluevl[2008ah : 

(false alarm probability) w We^^y/z, (13) 

where z is the height the maximum periodogram peak, and W = fmmJen is the normalized 
frequency bandwidth, with /max being the maximum frequency being scanned and T^ff being 
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the effective time series span (which is usually close to the actual time base). The applica- 
bility of this formula to the likelihood ratio periodograms was also discussed in ( iBaluevl . 
HoOSc). None of the peaks on the periodograms in Fig.[TO]possesses the false alarm proba- 
bility estimation less than ~ 20%. This means that no extra detectable periodicity is present 
in the RV data being used in our work. The only suspiciuous peak is shown by the peri- 
odogram of the residuals to the solution FA (the upper-right panel in Fig. \T0[ . This peak is 
close to the period of 26.6 days and possesses the false alarm probability of ~ 2%, if the 
latter is calculated for the period range P > 10 days (instead of P > 1 days). However, this 
peak is not present in other periodograms and thus seems to be a noise artefact exacerbated 
by inaccuracies of the model of the RV curve for this orbital solution and by aliasing (caused 
by uneven spacings of the RV data), rather than to reflect an RV oscillation induced by an 
extra hypothetical planet. 

There is also possibility that the putative additional planets are also trapped in a MMR 
with one of the two outer planets. For example, the RV oscillation of the hypothetical fourth 
planet having orbital period of Pd/2 w 450 days or Pc/3 w 650 days would be extremely 
difficult to extract from the synthetic RV curve. Such RV oscillation could be almost equally 
treated as a Fourier overtone associated with Keplerian oscillations induced by the outer 
planets (thus resulting in some change of their best-fitting eccentricities). Since the cur- 
rently available RV data for HD37124 do not allow reliable determination of the orbits in 
the system even for the three-planet configuration, the more complicated (and worse deter- 
mined) four-planet configurations were not considered here. It is worth emphasizing that 
there is no necessity to call for four-planet configurations of HD37124. All present difficul- 
ties connected with this system can be explained from the positions of data analysis, that is 
by the lack of the RV data, which is not sufficient to constrain reliably the non-linear triple- 
planet RV mod el having large number of degrees of freedom. For example, it was shown by 
iBaluevI ( l2008bl) that the unrealistically large formal estimations of the eccentricities can be 
interpreted as a result of their statistical (systematic) biasing. 

9 Conclusions 

In the paper, the full set of high-precision RV data available for the planetary system of 
HD37124 is analysed. The analysis involves different RV models and accounts for the re- 
quirement of the dynamical stability of the planetary system. The most likely orbital config- 
urations of the system, found in the paper, are split in four classes: 

1. Two outer planets 'c','d' are captured in the 2/1 mean-motion resonance and move on 
orbits with low or moderate eccentricities (e < 0. 15). The planets are far from an apsidal 
corotation resonance, but the whole system is still stable due to relatively low eccentric- 
ities. 

2. The planets 'c','d' are in the 2/1 mean-motion resonance and move on significantly 
elliptic orbits with moderate or high eccentricities. The system is stable, because the 
planets are locked in (or librate around) an asymmetric apsidal corotation. 

3. The planets 'c','d' are in the 5/2 mean-motion resonance. They move on elliptic orbits 
with high or moderate eccentricities. The planetary orbits are intersecting (or close to 
an intersection), but the system is stable, because the planets are locked in (or librate 
around) a symmetric antialigned apsidal corotation. This branch of solutions was also 
groped by Gozdziewski et al. (2006) basing on the Keck data only. 

4. The planets 'c','d' are not necessarily trapped in a mean-motion resonance. They move 
on orbits with relatively low eccentricities, which make the whole system stable. How- 
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ever, these solutions show larger scatter of the RV residuals and thus are less likely then 
solutions from other branches. When the eccentricies are fixed at small values, the ratio 
of the best fitting orbital periods of the two outer planets appears within a few per cent 
of the 2/1 MMR. The values of the period ratio Pc/Pd exceeding 2.3 — 2.5 are unlikely. 
However, the configurations with Pc/Pd < 2.3 are quite possible, bearing in mind the 
example of the planets b and e in the system of /J Ara, which are close to, but likely 
not tr apped in, the 2/1 MMR with Pi,/P^ w 2.1 dOozdziewski et all I2OO7I : IShort et aU. 
l2008h . 

The branches (2) and (3) are of a special interest, because we do not know examples of 
planetary systems with similar orbital configurations. The possibility of suc h orbital config- 
urations makes the algorithms of optimal scheduling of RV observations ( lBaluevll2008di ; 
[pprd, 2008) extremely tempting to use for this star. 

It was also shown that the RV residuals do not contain any detectable extra periodicity 
which could provide clear evidence for a fourth planet in the system. However, for the solu- 
tions close to the 2/1 MMR, there are evidences for an extra variation in the RV residuals. 
This variation can be explained by putative annual errors in the Keck data or, alternatively, 
by extra linear RV drift which could be induced by a distant unseen companion of the star. 
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A Obtaining ACR fits 

Here we describe the procedure of obtaining the ACR fits in more details. The algorithm that we are about 
to describe may be useful also for other planetary systems. But before we proceed, we have to choose some 
coordinate system. Actually, in the case of HD37124, it may be checked that the offsets of resulting orbital 
parameters, referenced in different coordinate systems, would be quite negligible in the sense of the RV fit 
quality. However the type of the coordinate system should be stated for the purposes of long-term integrations 
of the planetary system: the coordinate system used in the integration should match the given orbital elements . 
The actual choice of the Jacobi coordinates was motivated here by the fact noted in ('L issauer & River3 . l2001l ; 
[Lee & Peale, 2003), that it is Jacobi coordinate system in which the osculating orbital elements are mostly 
close to those obtained using the kinematic (Keplerian) RV model, especially when the system contains 
hierarchical planet pairs (like pairs of planets b-c and b-d in the case of HD37124). This property of the 
Jacobi coordinates may be useful, for instance, during a transition from multi-Keplerian to N-body model 
of the RV, that was (and probably will be) needed for some resonant planetary systeins after accumulating a 
sufficiently long observation time span. 

In the Jacobi coordinates, the Hamiltonian of a system with N planets looks like 



^-E^-^ - E (14) 



l<i<j<JV 



where G is the gravitational constant, m, are the planetary masses, r,- are the astrocentric distances of the 
planets, r,j are the distances between the planets, pj are the Jacobi momenta andmj = "iiCEjCQnij) / {Y,'j=a"'^j) 
are the Jacobi masses. Now we need to split H4i in the Keplerian part that we consider as unperturbed one 
and in the part that we consider as perturbational function. This may be done non-uniquely. We adopt the 
splitting H = I «K.p., - R„ where 

[pf GM.^i \ „ ^ [ Mj mo M,A 

^Kep.- = 1 ^ - I : = 1 E — + — - — ) ■ (1^> 



j=i+l '•] 



Here rj is the length of the Jacobi radius- vector for the planet, M,_ i = I^^^'q r. 
mass and of the masses of all planets being interior with respect to the given planet. Therefore, //Kcp.i is 
chosen so that the con'esponding unperturbed Keplerian orbit is referenced to a fictitious central body having 
mass Mj-\. We will use only the second-order approximation of the Hamiltonian. In this approxiination, we 
represent S, = Y!j=i+i ^ij + &{m^), where 

Rij = Gm.mj ( - ^ - \] . (16) 

Here r' are the Jacobi p osition vectors. This approximation for the case of two planets inay be found in 
iGerassimov et al.Lll99d § 4.2), and the extension to A' > 3 planets is straightforward. 
Let us define the following function of six input parameters: 

^(a,einn,eout,.^,m,A„ut,4(B) = ^ - """"j ^— . (17) 

rinn r^yt I rjj^j ^out 

Here, the vectors T\^^ and rout describe positions of abstract 'inner' and 'outer' planets moving on Keplerian 
orbits. The Keplerian orbital elements marked by the subscript 'inn' refer to the orbital elements of the inner 
planet, and quantitties marked by the subscript 'out' correspond to the outer planet. The semi-major axis of 



27 



the inner abstract planet is set to unit, and the semi-major axis of the abstract outer planet is equal to a. The 
parameter A(0 = COout — <'Km\ represents the difference between the longitudes of periastra. Then the obvious 
identity Rij = "'""^ ^(a ; /a, , e, , ej , A, . Ay , (Oj — (O-,) holds true. 

As usually, the Keplerian parts of the Hamiltonian depend only on the Delaunay actions L, (and, as on 
parameters, on the masses of the planets) and remain unchanged by any averaging. Let us now restrict our 
attention to the motion of the (possibly) resonant planets HD37124 c and d, neglecting the influence of the 
innermost planet b. The perturbational function describing the interaction of the planets c and d is given 
by "a^"''' .^(ct = ac/aj,ej,ec,^d,^,Aco = cOc — coj). The corresponding averaged perturbational function 

is then given by {^){oc,ej,ec,sj, .sv ) , where the averaged function {M) depends on only five input 

parameters a , , eom , -^inn , ^om : 

f^" d9 

JO 271 

Note that the last term in the expression I17t is averaged to I /a. This term reflects only the fact that the 
osculating Keplerian orbits are refered to different fictitious central masses. In fact, we need to average only 
the resting classical expression of the perturbational function: 

{.'%) = (. i (19) 



' out 



a 



We perform the averaging U8|19) by means of numerical integration tools, as it was proposed bv lMichtchenko et alj 
12006). This way of numerical averaging of the Hamiltonian is very easy to implement and simultaneously is 
quite rapid and precise. In the same way, we can calculate various derivatives of the averaged function (^), 
that we will need below. For example, 

^ d{R) ^ / rinn ■ rout - r^^, ^ ^ ''inn - Tout \ _|_ _L ^20) 

da \ Irinn-routh"* r'^, / a' 

where the averaging can be again performed numerically. 

In the next step, we need to build in our RV fitting algorithm the four equality bounds dlHcd) /d{IcJd) = 

and (9(//(.rf)/(9(.sv,.V(/) = d{R^j) /d{sc,Sd) = 0, which determine the location of the ACR. The second 
pair of equations reflects the fact that sj and should correspond to the extremum value of {M){a = 
ac/aj,ej,ec,sj,Sc) (considering the values of e^d and a fixed). This condition can be used to construct the 
angles s^.d as functions of the eccentricities e^- d- Sd = s*^^{a,ed,ec) and .sv = a, e^, fie), where the func- 
ti™** *,*nn("''^inn,eout) and C («, gjnn , fiou, ) provide an extt-emum to (.*>(a,ei„n,eou,,ii„n,ioi,t) given fixed 
0!,eitini£out- When the ACR is symmetric, the functions s*^^ and can be found easily (actually, they ap- 
pear to be constant in this case). However, for an asymmetric ACR, we have to locate the values of .v^ and 
Sd numerically. Eventually, this numerical procedure is still sufficiently rapid and can be implemented as a 
'black-box' subroutine, which returns the values of s*^^^ and .s*^ for input values of einn, eout, and ct- Note that 
since we consider here only MMR solutions, a ^ ao = (p/qf'l^ with a eiTor of &(mcjlMf,). When calcu- 
lating the resonant angles, we can quite neglect such en'ors and put simply a = Oq. Such errors in resonant 
angles will not produce significant changes in the dynamics of the planetary system. The resulting values of 
Sc and Sd and the definitions jlOt can be used in the work of the fitting algorithm to express tOc and Wd via 
the other free parameters: , , , . 

Another pair of constraints, d{Hcd) I dijc-jd) = 0, can be transformed to a more convenient (and equiv- 
alent) form, involving partial derivatives over the Delaunay actions L^ d and Gc.d = L^ dTlcd (here rj^j = 

1 — ^). The first transformed equation, d {Hcd) / dGc = d (H^d) /dGd, reflects the coincidence of the secular 
drifts of the orbital periastra. After neglecting insignificant errors of the order of the planetary masses, this 
equation can be simplified to 



^ = aop(ao,erf,ec) with p{a,e,„„,e„„t) = ji'fwj'^'"" > (21) 

f^c '5(.«>/'5i?out (.i„„,.ou.) 

■^out=SoLi,(Pinn''out> 

where we note that ^ = — i ^ . The equality )2U can be used to express the RV semi-amplitude of one of 
the resonat planets via the RV semi-amplitude of another one and via the orbital eccentricities. The second 
transformed equation, pd(Hcd) /dLc = qd(Hcd) jdLd, reflects the vanishing of the secular drift of the critical 
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angle pic — qld (where / are the planetaiy mean anomahes). This should provide the long-term constancy of 
the planetary conjunction positions. After some simplifications, this equation may be rewritten as 



v(a,ein„,eoiit) - 



P 

q 



■-21 a 



Pd_P 
Pc q 



1 '"c / \ 



with 



da 



1 + 

q 



(22) 



This equation introduces a error of the second order only, that is ff ((m/M,t)^). Actually, specifically to the 
system of HD37124, the small ff(m/Mi) deviation of the period ratio from the exact resonance does not affect 
significantly the quality of the RV fit (moreover, the values of V for the ACR fits from Table [5] appeared less 
than 0.1). Almost the same value of the RV r.m.s. could be obtained for the simplified equation Pc = ^Pj- 
However, it is this 0(m/Mi,) period ratio deviation that determines the secular drift of the critical angle 
pic — qld and the ACR state of the system. Therefore, in general case it is necessary to take this deviation 
into account when constructing the ACR fits. The equation j22) can be used to express the osculating orbital 
period Pc or Pj via the resting free variables. Note that the deviation of the ratio of the osculating semi-major 
axes can be written down (to within the first order) as 



a _^ Mc 2v{ao,eci,ec) + l 
Oo M* 3 ' 



(23) 



Therefore, for any given values of the parameters ec,ed,Xc,Xd,Kci,Pd we can obtain the ACR val- 
ues of the parameters (OcCOd-Kc with a error of ff{mc^d/Mi,) and the ACR value of Pc with a en'or of 
6" {{I'nc.d If^t-kY') ■ This allows us to obtain the best-fitting orbital solution, which is sufficiently close to an 
ACR state. 



